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A  MATHEMATICAL  METHOD  FOR  INVERSION  IN 
ATMOSPHERIC  REMOTE  SENSING 

INTRODUCTION 

The  importance  of  Fredholm  integral  equations  of  the  first  kind  is 
well  known  to  applied  mathematicians.  One  reason  for  this  is  that  they 
occur  in  almost  every  area  of  the  exact  sciences  either  in  the  formulation 
and  solution  of  theoretical  problems  or  in  the  decomposition  of  sets  of 
measurements  to  obtain  desirable  values  and  other  types  of  information. 
Examples  of  their  occurrence  include  convolution  equations,  integral 
transforms,  instrumental  broadening,  and  inversion  of  indirect  measure¬ 
ments,  in  addition  to  convenience  in  expressing  some  physical  laws.  In 
areas  of  remote  sensing  amenable  to  quantitative  description,  the 
opportunities  for  their  application  are  particularly  numerous.  Indeed, 
for  successful  solution  of  many  remote  sensing  problems  they  are  essential 
and  unavoidable;  yet  they  remain  unused  because  of  numerical,  analytical, 
and  other  difficulties,  and  because  researchers  in  this  area  often  have 
backgrounds  in  inexact  and  descriptive  sciences  and  neither  understand  nor 
appreciate  the  importance  of  integral  equations. 

This  has  not  been  true  in  the  area  of  atmospheric  inversion,  where  as 
early  as  1959  Kaplan  suggested  that  It  was  possible  to  use  remote  Infrared 
measurements  to  infer  the  structure  of  the  atmosphere.  In  spite  of 
considerable  skepticism,  which  was  based  partly  on  well-founded  physical, 
experimental,  and  mathematical  difficulties,  progress  has  been  substantial 
and  a  limited  degree  of  success  has  been  achieved. 

From  an  analytical  point  of  view,  the  difficulty  in  solving  Fredholm 
equations  of  the  first  kind  often  arises  because  the  equation  transforms 
or  maps  the  manifold  on  which  the  unknown  functions  are  defined  onto  a 
different,  smaller,  or  more  restrictive  manifold.  In  such  cases  there  may 
no  longer  be  a  unique  solution  because  the  mapping  of  one  function  space 
onto  the  other  is  "many  to  one."  Well-known  examples  include  band-limited 
convolution  products  and  integral  transformations  of  discontinuous 
functions  using  continuous  kernels. 

For  problems  where  a  numerical  (not  an  analytical)  type  of  solution  is 
needed,  the  Integral  equation  formulation  constitutes  an  Incorrect  or  ill- 
posed  problem  in  the  sense  of  Hadamard,  i.e.  a  unique  solution  depending 
continuously  on  the  transformed  or  unknown  function  does  not  exist. 
Nevertheless,  by  redefining  the  problem,  correctness  can  be  established  in 
other  ways  and  solutions  can  be  constructed.  In  addition  to  the  problem 
of  correctness,  severe  limitations  in  the  amount  of  available  data,  errors 
in  data,  and  limitations  in  the  numerical  methods  used,  combined  with 
usual  computational  inaccuracies  and  possible  instability,  have  led  to 
solutions  that  were  meaningless  and  even  physically  impossible  for  many 
methods  that  have  been  tried. 
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The  method  to  be  proposed  here  attempts  to  overcome  these  difficulties 
by  introducing  certain  types  of  discretization  and  by  restricting  the 
function  space  of  the  solution.  Only  the  formal  aspects  are  presented  and 
the  method  is  general  since  it  leaves  many  choices  to  the  ingenuity  of  a 
user.  It  is  related  to  some  other  methods  that  have  been  used,  and  the 
discussion  should  help  to  elucidate  some  of  the  problems  encountered  in 
using  them. 


THE  TRANSFER  EQUATION 

Before  proceeding  to  the  mathematical  method  for  its  solution,  a 
discussion  of  the  transfer  equation  will  be  given.  Although  this 
formulation  is  well  known,  the  brief  discussion  will  define  the  problem  in 
physical  terms  and  show  how  the  transfer  equation  is  expressed  as  a 
Fredholm  equation. 

It  is  assumed  that  the  transfer  of  energy  within  the  atmosphere  occurs 
by  radiation  through  a  clear  column  of  air,  that  clouds  are  not  present, 
and  that  scattering  processes  are  unimportant.  Furthermore,  a  local 
condition  of  radiative  thermodynamic  equilibrium  is  assumed  so  that  the 
ratio  of  the  coefficient  of  emission  to  that  of  absorption  is  given  by  the 
Kirchhoff -Planck  equation 


B( v,  T)  - 


2hv3 

r»2 


where  v  is  the  frequency  of  the  radiation  and  T  is  the  absolute 
temperature.  This  represents  a  convenient  abstraction  for  situations 
where  continuous  spectra  are  encountered  and  where  no  correlation  exists 
between  the  particular  emissions  and  absorptions  that  participate  in  the 
multiplicity  of  processes  Involved. 

If  p  is  used  for  pressure  level,  the  temperature  profile  is  defined  by 
T(p)  corresponding  with  measured  radiance  Intensity  I(v),  and  an  initial 
approximation  is  assumed  to  be  given  by  T(p),  I(v).  The  quantity 


A  l(v)  -  l(v)  -  I(v) 


is  then  expressable  by  a  different  form  of  the  radiative  transfer  equation 


&I(v)  -  «B(s)  Ty(s)  +  /  AB(p)  dTy(p) 

Tv(s) 


where  Tv(p)  is  the  transmittance  of  the  atmosphere  from  pressure  level  p 
to  the  top  and  where  s  refers  to  the  surface  pressure  level.  Furthermore 
SB  and  AB  are  given  by 


6B(s)  -  B(vr,  T(p))  -  B(vr,  T(s)  )  (4) 

AB(s)  -  B(vr,  T(p) )  -  B(vr,  T(p) )  (5) 


where  T(p)  is  the  assumed  surface  temperature  and  has  been  a  particular 
source  of  difficulty.  By  noting  that 


J1  dTv(p)  -  1  -Tv(s)  (6) 

V«) 

both  terms  of  the  transfer  equation  (3)  can  be  combined. 

AI(v)  -  f1  iAB(p)  +  6B(a)  }dtv(p)  (7) 

*v(s)  l-tv(s) 


It  la  to  be  expected  that  a  correspondence  (not  necessarily  one  to  one) 
exists  between  the  bracketed  term  and  some  function  of  the  pressure 
profile 


h(p) 


tv(s) 

1  ”  Tv(s) 


(8) 


so  that  the  transfer  equation  assumes  the  customary  form  for  this  problem 


where 


dtu(p) 

AI(v)  -  -J  x(p)  -  dp  (9) 

o  dp 

x(p)  -  AB(p)  +  6B(s)  h(p)  (10) 


In  this  form,  the  radiative  transfer  equation  (9)  is  a  Fredholm  integral 
equation  of  the  first  kind  to  be  solved  for  x(p)  where  the  boundary 
term  6b(s)  h(p)  is  to  be  subtracted  to  obtain  AB(p).  Since  the  function 
h(p)  is  not  known,  it  could  be  constructed  by  integrating  the  transfer 
equation  for  many  known  widely  divergent  profiles  at  several  wavelengths 
and  using 


h(p)  *  [x(p)  -  AB(p)]/fiB(s) 


(11) 


but  it  more  often  must  be  obtained  by  using  surface  transmittances . 
Different  functions  h(p)  are  characteristic  of  atmospheres  of  different 
temperature  and  humidity  type. 


THE  MATHEMATICAL  INVERSION  METHOD 

In  more  conventional  mathematical  notation,  the  equation  of  transfer 
(9)  can  be  written  in  the  form 

b 

f(x)  -  /  K(x,  «  u(5)  dC  (12) 

a 

where  the  function  f  is  given  by  a  series  of  measurements  for  fixed  values 
of  x  and  where  the  kernel  K  is  considered  to  be  known.  The  problem  of 
inversion  consists  of  recovering  the  unknown  function  u  to  within  some 
reasonable  limits.  In  atmospheric  inversion  problems,  the  kernel  function 
is  smooth  and  limits  the  ability  to  recover  higher  frequency  components 
resulting  in  an  ambiguity  (or  in  some  methods,  a  so-called  "inherent 
instability")  in  the  inversion  for  f. 

In  order  to  effect  an  inversion,  a  set  of  basis  functions  ut(x)  is 
selected  to  span  a  function  space  of  low  dimension,  but  capable  of  giving 
a  representation  of  the  solution  u(x)  in  the  form 

N 

u(x)  -  l  03) 

i-1  1 

This  is  suggested  by  the  linearity  of  the  integral  equation  (12),  but  the 
basis  functions  need  not  physically  represent  actual  atmospheric  profiles. 
They  need  only  be  capable  of  being  summed  to  give  a  good  approximation  to 
the  profile  being  sought.  The  solution  will  consist  of  determining  the 
constants  a^. 

From  the  basis  functions  uj,  a  new  set  of  N  functions  f^(x)  is  to  be 
computed  by  performing  the  quadratures  indicated  by 

b 

fj(x)  •  /  K(x,  5)  u1(5)  d£  (14) 

a 
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for  the  values  of  x  at  which  measurements  have  been  made. 

Another  set  of  N  functions  $j(x)  Is  also  to  be  constructed,  at  the 
same  values  of  the  argument  x,  as  linear  combinations  of  the 
functions  f^(x) 


N 

4>j(x)  -  l  fk(x)  (15) 

where  the  constants  ctj^  are  determined  by  imposing  the  condition 

b 

/  <frj(x)  f£(x)  dx  -  (16) 


where  the  Integral  implies  a  numerical  quadrature  over  the  appropriate 
values  of  x.  Multiplying  both  sides  of  equation  (15)  by  f&(x)  and 
performing  a  similar  quadrature  leads  to 

b  N  b 

/  <|>  (x)  f£(x)  dx  -  l  /  ffe(x)  f£(x)  dx  (17) 

a  J  k-1  J  a 

which  becomes 


(18) 


when  equation  (16)  is  introduced  and  the  constants  Pki  are  defined  by 
equation  (19). 


\i  - 1  fk(x)  fiM  d* 


(Id) 


The  quantities  8^4  are  compatible  by  the  numerical  quadrature  indicated  In 
equation  (19)  since  the  functions  ffc(x)  are  known.  Then,  ajk  satisfying 
equation  (18)  can  be  found  by  taking  the  matrix  Ojk  to  be  the  inverse 
of  {djfc}  where  it  may  be  convenient  to  make  use  or  the  fact  that  both 
matrices  must  be  symmetrical.  The  functions  +j(x)  can  now  be  obtained  by 
use  of  equation  (15).  It  will  be  found  that  they  are  not  needed  for  a 
final  computation,  but  are  useful  In  the  derivation  of  the  solution. 
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Inserting  equation  (13)  into  the  integral  equation  (12), 

b  N 

f(x)  =  /  K(x,  5)  l  a  u  (C)  dC 


N  b 

f(x)  =  la/  K(x,  C)  u  (5)  dC 


ra-1  a 


it  is  found  by  using  equation  (14)  that 


f(x)  -  l  a  f  (x). 
m-1  m  m 


Multiplying  both  sides  of  equation  (22)  by  4n(x)  and  integrating 
(performing  a  quadrature)  yields  a  method  for  computing  the 
coefficients  an. 


b  N  b 

/  $n(x)  f(x)  dx  *  l  a^  /  <frn(x)  fm(x)  dx  (23) 

a  m-1  m  a 


£  a  6 
,  m  nm 
m-1 


/  $n(*)  f(x)  dx  -  aR 

a 


The  insertion  of  equation  (15)  for  $n(x)  into  equation  (25)  leads  to  the 
proper  expression  for  computing  the  coefficients  aj  in  equation  (13), 

b  N 

a  -  /  I  a  f  (x)  f(x)  dx  (26) 

n  a  k-1  nk  k 


V  !  Vx)  f(x)  dx 

k-1  a 
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where  Y^  is  defined  as  follows 


b 

Yk  -  /  fk(x)  f(x)  dx  (29) 


The  final  solution  can  now  be  written  by  replacing  aj  in  equation  (11) 
with  equation  (28). 


u<x>  '  £,  J,  V  Tk  (30> 

n-1  k-1 


Recapitulating,  the  steps  in  finding  the  solution  consist  of  performing 
the  quadratures  indicated  in  equations  (14),  (19),  and  (29)  and  inverting 
to  obtain 

It  should  be  remarked  that  the  same  type  of  quadrature  formula  can  be 
used  for  most  cases  involved.  Also,  N  must  not  exceed  the  number  of 
points  at  which  the  function  f(x)  is  given  if  the  matrix  t8nji}  is  to  be 
successfully  inverted. 


CONCLUSION 

A  method  of  solving  Fredholm  integral  equations  of  the  first  kind  for 
use  in  the  numerical  inversion  of  remotely  sensed  radiance  measurements 
has  been  presented.  Use  of  this  method  requires  the  specification  of  a 
basic  set  of  functions  from  which  the  solution  can  be  written  as  a  linear 
superposition.  It  does  not  impose  further  unnecessary  conditions,  such  as 
the  orthogonality  often  used  in  obtaining  such  solutions.  This  should 
allow  a  good  representation  among  the  physically  realizable  solutions  from 
a  relatively  small  number  of  basis  functions.  This  is  important  since  the 
information  available  for  the  function  to  he  inverted  is  usually  severely 
limited  by  physical  restrictions  on  the  number  and  types  of  measurements. 

The  general  mathematical  aspects  of  the  problem  and  a  brief  derivation 
of  the  mathematical  equation  from  physical  considerations  have  been 
presented. 
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